
#include <math.h>
#include <process.h>
#include <malloc.h>
#include "rebound.h"

double e_init;
double tmax;
float* _values;
int _length;
int _traceCount;
int index;

__declspec(dllexport) void Test(
	float* values,
	int length,
	int traceCount
	)
{

	struct reb_simulation* r = reb_simulation_create();

	// Setup constants
	//r->dt = 0.01;             // in days
	r->G = 1;// 1.4880826e-34; // in AU^3 / kg / day^2.
	//# All units entered in SI(kg, m, s)		G_SI = 6.67408e-11
	r->ri_ias15.min_dt = 0;
	r->ri_ias15.max_dt = 1;
	//r->ri_whfast.safe_mode = 0;             // Turn off safe mode. Need to call reb_simulation_synchronize() before outputs. 
	//r->ri_whfast.corrector = 11;            // 11th order symplectic corrector
	//r->integrator = REB_INTEGRATOR_WHFAST;
	//r->collision = REB_COLLISION_DIRECT;
	//r->heartbeat = heartbeat;
	r->exact_finish_time = 0;     // Finish exactly at tmax in reb_simulation_integrate(). Default is already 1.
								  //r->integrator         = REB_INTEGRATOR_IAS15;        // Alternative non-symplectic integrator
								  // Initial conditions

	for (int i = 0; i < length; i++) {
		struct reb_particle p;
		p.x = values[i * 8 + 0];  p.y = values[i * 8 + 1];  p.z = values[i * 8 + 2];
		p.vx = values[i * 8 + 3]; p.vy = values[i * 8 * 8 + 4]; p.vz = values[i * 8 + 5];
		p.m = values[i * 8 + 6];
		p.r = values[i * 8 + 7];
		p.ap = (void*)i;
		reb_simulation_add(r, p);
	}
	_values = values;
	_length = length;
	_traceCount = traceCount;
	index = 0;
	//reb_simulation_integrate(r, tmax);

	int index;
	for (index = 0; index < traceCount; index++) {
		int itr = 0;
		while (r->t < index*1.0) {
			reb_simulation_step(r);
			itr++;

			int i;
			int j;
			struct reb_particle* p1;
			struct reb_particle* p2;
			int found = 0;
			for (i = 0; i < r->N; i++) {
				for (j = 0; j < r->N; j++) {
					if (i != j) {
						p1 = &(r->particles[i]);
						p2 = &(r->particles[j]);
						if (p1->m != 0 && p2->m != 0) {
							double dx = p1->x - p2->x;
							double dy = p1->y - p2->y;
							double dz = p1->z - p2->z;
							double dist = sqrt(dx * dx + dy * dy + dz * dz);
							double dist2 = pow(p1->m, 0.33333333333) + pow(p2->m, 0.33333333333);
							if (dist < dist2) {
								found = 1;
								break;
							}
						}
					}
				}
				if (found) {
					break;
				}
			}
			if (found) {
				if (p1->m > p2->m) {
					p1->m = p1->m + p2->m;
					reb_simulation_remove_particle(r, j, 0);
				}
				else {
					p2->m = p1->m + p2->m;
					reb_simulation_remove_particle(r, i, 0);
				}
			}
			/*
			if (itr > 100000) {
				//r->particles[i].

				for (int i = 0; i < _length; i++) {
					printf("p:%d x:%f y:%f z:%f \r\n", i, r->particles[i].x, r->particles[i].y, r->particles[i].z, r->particles[i].m);
				}
				return;
			}
			*/
		}
		
		for (int i = 0; i < r->N; i++) {
			int idx = (int)r->particles[i].ap;
			_values[_length * 8 + index * _length * 4 + idx * 4 + 0] = r->particles[i].x;
			_values[_length * 8 + index * _length * 4 + idx * 4 + 1] = r->particles[i].y;
			_values[_length * 8 + index * _length * 4 + idx * 4 + 2] = r->particles[i].z;
			_values[_length * 8 + index * _length * 4 + idx * 4 + 3] = r->particles[i].m;
		}
		printf("step idx:%d x:%f t:%f dt:%f itr:%d\r\n", index, r->particles[0].x, r->t, r->dt, itr);
		//r->dt = 0.1;
	}
}





/*

for cpp
extern "C" __declspec(dllexport) void Test(
		float* posx,
		float* posy,
		float* posz,
		int length
	)
{
	posz[0] = posx[0] + posy[0];

	struct reb_simulation* r = reb_simulation_create();


}
*/